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SHEAR  STABILIZATION  OF  DRIFT  WAVES  IN  NONCIRCULAR 
CROSS  SECTION  AXISYMMETRIC  CONFIGURATIONS 

I.  INTRODUCTION 

The  role  of  shear  in  the  stabilization  of  drift  waves  has  received 
a great  deal  of  attention  * (For  example,  the  treatment  of  the 

q 

dissipative  trapped  electron  mode  instability  has  recently  been 
studied  because  of  its  interest  for  tokamak  transport.)  In  the  most 
simplified  approaches  the  growth  rate  may  be  expressed  as 

Y = Y£-Ys 

where  Yp  is  a mode  growth  term  usually  due  to  an  electron  dissipation 
mechanism  (e.g.,  scattering  of  trapped  electrons  for  the  dissipative 
trapped  electron  mode,  electron  Landau  resonance  for  the  universal 
instability^,  finite  electron  thermal  conductivity  for  the  drift 
dissipative  mode),  and  vs  is  a shear  induced  damping  rate.  Shear 
stabilization  results  when  vs>  y,  . The  calculation  of  y was  given  by 
Pearlstein  and  Berk^  for  slab  geometry  [b  = B (zQ  + yo  x/l  ) where 
L$  is  the  "shear  length")  in  their  study  of  the  universal  instability. 
They  took  the  electrostatic  potential  in  the  form  f(x)  exp(ik  y)  and  found 
a mode  localized  about  x=0,  which  propagates  energy  outward  in  x until, 
as  a result  of  the  shear,  ^/v.  %|k„|  - k (x ( / L s , where  the  wave  energy 
is  damped  via  ion  Landau  resonance  ( v . the  ion  thermal  velocity).  This 
calculation  is  trivially  extended  to  a circular  plasma  <.  nlumn  equilibrium 
which  models  a large  aspect  ratio  tokamak  of  rum  i , ! if  we  impose 

mode  periodicity;  along  the  column  length,  z - ft  t / Th.  mason 

for  the  essential  similarity  between  the  stub  md  . i-.  . 
geometries  is  that  they  both  have  two  direction,  of  un  • Uv,  v .• 

Note:  Manuscript  submitted  May  2,  1977. 
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in  the  slab  model,  and  0 and  z in  the  circular  cross  section  case.  On 

the  other  hand,  there  is  much  interest  in  tokamaks  with  vertically 

elongated  cross  sections  due  to  their  potential  advantage  with  respect 

1 ? 

to  achieving  higher  3 and  higher  current  densities  . In  addition, 
multipole  devices  also  exhibit  magnetic  surfaces  with  noncircular  cross 
section.  In  these  devices  the  geometrical  0 symmetry  of  circular  cross 
section  is  lost. 

It  is  the  purpose  of  this  paper  to  obtain  the  shear  induced 
damping  rate  of  drift  waves  in  such  noncircular  cross  section  devices. 

We  make  the  following  simplifying  assumptions:  (i)  large  aspect  ratio, 

(ii)  electrostatic  perturbation,  (iii)  modal  wavelengths  larger  than 
Pi(pi  is  the  ion  Larmor  radius),  and  (iv)  a longitudinal  magnetic  field 
much  larger  than  the  poloidal  field. 

Section  II  describes  the  geometry  and  obtains  the  fundamental 
two  dimensional  partial  differential  equation  which  governs  this  problem. 
Section  III  presents  the  solution  of  the  fundamental  equation  for  the 
shear  induced  damping  under  the  assumption  that  mode  rational  surfaces 
are  decoupled  and  discusses  this  assumption.  As  an  example,  these 
results  are  applied  to  a particular  tokamak  equilibrium  in  Section  IV, 
and  it  is  shown  that  the  shear  induced  damping  rate  is  very  insensitive 
to  the  amount  of  vertical  elongation  (at  least  for  the  example  given). 


■wnn mi 


Related  work  has  been  reported  by  Connor  and  Hastie^  who  use 
a different  method  and  do  not  include  plasma  current  (applicable  to 
levitrons  and  multipoles  but  not  tokamaks).  The  main  difference  with 
their  result  is  that  the  da/d'Jj  terms  in  the  frequency  shift  d,  given  in 
our  Eq.  (14),  would  be  absent  since  without  plasma  current  a is  constant. 

O 

Talor  has  treated  a model  equation  to  determine  the  effect  of  strong 
coupling  between  mode  rational  surfaces.  Here  we  only  treat  situations 
where  mode  rational  surfaces  are  decoupled  by  ion  Landau  damping 
(cf.  discussion  at  the  end  of  Section  III). 
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II. 


GEOMETRY  AND  FUNDAMENTAL  EQUATION 


We  assume  a noncircular  cylinder  of  length  2ttR  with  periodic 
boundary  conditions  applied  at  the  column  ends  (cf.  Fig.  1(a))  to 
simulate  a large  aspect  ratio  axisymmetric  device.  The  coordinates  are 
( ip,x,z ) and  are  illustrated  in  Fig.  1,  where  constant  ip  corresponds  to 
a magnetic  surface  d ip = B^dl^  , x specifies  the  crossectional  position 
of  a point  on  a flux  surface  dx=otB  dl  ,z  is  the  length  along  the  cylinder, 
B^  is  the  poloidal  component  of  B (component  X to  z),  and  dl^  and  dlx 
are  differential  lengths  in  the  <p  and  x directions.  The  function  a is 
chosen  to  make  the  (ip,x)  coordinates  orthogonal, 

Vx«  B =0.  (1) 

In  order  to  see  that  this  is  so  assume  ip  and  x are  orthogonal  as  in 
Figure  1 (b)  and  take  the  line  integral  § a S^dl  around  the  closed  path 
(a,b,c,d)  indicated  in  Figure  1(b).  We  see  that  for  an  orthogonal 
coordinate  system  /aB^-dl =0  and  so  (1)  follows. 

We  now  set  down  the  basic  equations  describing  the  system.  In 


doing  this  we  make  the  following  assumptions  to  simplify  the 
(i)  u-»k„vi,  (ii)  Tg»T i , (iii)  w«k„ve,  (iv)  k±2p.2«l. 

analysis 

"i  = "e  = 

(2) 

ne  = N(iJ;)  e$/Te, 

(3) 

-iton.j+  V-  (n^v)  = 0, 

(4) 

v = vx  + v„ 

(5) 

v±  = - B“2V<t>  XB  + iu)(Bw,.  )_1vx4> 

(6) 

4 


-■MV,, 


Cl 


B • V4> 


(7) 


where  = eB/M^  and  E = - VO. 


Eliminating  h and  v and  using  tokamak  ordering  ( B B ) < < 3 and  taking 
0 'v  exp(ikz),  k=n/R,  we  obtain  a single  second  order  partial 

differential  equation  forOin  terms  of  the  coordinates  ip  andx  . The 
details  leading  to  this  equation  appear  in  Appendix  A.  The  result  is 

c 2 

9 . s , _ n 2 3 


icT_  o „ c.  n „ 2 


eB 


A* 


cTp  w 

"eB  aB-L  ( ; 


“B 

2/313 


"L  V 3(1/  a 3ip  3X  01  3x  )]  * °’  ^ 

1/2 

where  ^=eO/Te  and  c =(Tg/M.)  ' . Equation  (8)  is  our  basic  equation 
to  be  studied  in  detail  in  the  next  section.  The  first  t' rm  in  (8)  is 
due  to  ng;  the  second  term  is  due  to  the  ion  density  perturbation  produced 
by  the  E x B drift  and  background  density  gradient;  the  third  term  is 
a result  of  finite  ion  inertia  along  field  lines;  and  the  last  two  terms 
in  (8)  represent  finite  ion  inertia  across  B (the  polarization  drift). 

If  the  destabilizing  sources  of  electron  dissipation  were  included 


Eq.  (3)  would  take  on  the  general  form  he=N(\J/)  |l+iivi)  (e<t/Tc)  where  M 

A 

can  be  an  operator.  Different  forms  of  M result  from  different  electron 

A 

dissipation  mechanisms  (as  discussed  in  Section  I).  Here  we  set  M :D 
(so  that  no  instability  mechanism  is  present)  and  solve  for  the  mode 
structure  and  shear  induced  damping  rate.  The  effect  of  M for  specific 
instabil ities  can  then  be  included  as  a perturbation  on  this  basic 
solution  as  demonstrated  in  References  4 and  5. 
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III. 


ANALYSIS  OF  BASIC  EQUATION 


To  analyze  Eq.  (8)  we  expand  v in  a set  of  orthogonal  basis 
functions  um(\,ij,) 


* (x.'l')  =T.  u (x,t^)  ^ (4,) . 

mm  m 


(9) 


It  is  possible  to  choose  the  un]  so  that  the  first  three  terms  in  (8), 
which  are  larger  than  the  last  two,  are  diagonal  for  all  values  of  ij;. 
This  will  be  the  most  convenient  choice  for  the  um.  We  choose  (nict  B^) 

um  - ex>>  [-1  „?STkB  !'  n'1  “*]  • »°) 


so  that 


$ u un  dx/n  = 6 tfdx/n, 
m p pirr 


(11) 


where  6 m is  the  Kronecker  delta,  * denotes  the  complex  conjugate,  and 
the  safety  factor  is  q ( ip ) = (2™)  ^kB^ii'^dx.  Note  that  the  um(x,^)  also 

i 

satisfy  the  periodicity  condition  um(x ,’p)=um(x+jFdx  slp)  on  all  magnetic 
flux  surfaces.  If  we  represent  Eq.  (8)  symbolically  as  Lv>  =0  then  use 
of  (9)  converts  (8)  to  a matrix  operator  equation 

r ■ o,  iio) 

where  each  of  the  Lpm  is  a second  order  ordinary  differential  operators 
in  p gotten  from  the  relation 


V'»  ■ I**  (“6lr‘  UP  L(um-',m>l  Ud'(“Bi)"ir1- 


(ID 


From  (11)  and  (8)  we  obtain 


= D m + 

pm 


N m + P m 
pm  pm 


(12) 


where 


eB 


|qW"m/nl  }’ 

(ijw)  Q(n)  I*”1'2  exp  1(Fp  - Fr)d,il 


P =-— p—  — — [£n  ^ d x ] fa  ^ exp  i (F  - F ) j 

pm  eB  wci  F 1 p m j d7 


1 3a 

t 3Fm] 

d f/81 

3Fm\)  j 

a 9ty 

+ 2 ‘ 

W j 

*T  Lv 

3*  j 3*  Jj  j 

where  Fm  = kBm/(nq(iJ>) ) /'n-1  dx, 

' w = -ckT  e-19«.n  N/9tK  and  D . N „ and  P „ represent  respectively 
*e  e pm  pm  pm 

the  first  three  terms  of  (8),  the  fourth  term  in  (8),  and  the  last 
term  in  (8).  Note  that  Dpm  is  diagonal  while  Npm  and  Ppm  are  non- 
diagonal, and  Ppm  is  a second  order  differential  operator  in  ijj  . Here 
we  only  consider  the  case  where  it  is  valid  to  consider  the  matrix 
operator  as  diagonal.  We  discuss  the  conditions  under  which  this  is 
valid  later.  Ignoring  the  off  diagonal  components  of  we  have  Lmrn*m=0 
which  using  the  approximation  q (i|j ) n m/n  + q'(^0)  (<l<  - <k0) 

(q’  = dq/dijj  and  q (»po ) = m/n)  may  be  written  as 
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i 

< 


) <i2  -ci  2'1 

) d7  ■/  ' ‘ 


1 - — d 

a) 


-( 


“•V  % « - v'][v» 

^ V0 


(13) 


where  d is  a frequency  shift  term  due  to  the  finite  cross-field  ion  inertia 
terms 


. s 1,2  „2  /n-2v 

d 9 } k B <B±  > + <Bx> 


-1 


Cl 


,3F  ,2 


<Bi  (w)  > <Bi> 


,9F_.V  2 


- <»i  (ri”)>‘  - 4 K ddf  >)  + ? <Bi> <*B1  ~r>l j , <14> 


the  flux  surface  averaqe  of  a function  6 is 
<G>  = 


ffc  B^dl^/fB^1  dlx 


J 


and  the  new  dependent  variable  has  been  introduced  to  eliminate  the 


first  derivative  term  in  P , 


*m  = exP  ( -/  9d^)v‘m(v)  . 

3F 

9 ( 4’)  = (/a_1dx)_1  (-*  j jfa_1dx  + ifa"1  “■  dx). 

Equation  (13)  is  analogous  to  the  equation  of  Pearlstein  and  Berk* . 
Looking  for  solutions  satisfying  an  outgoing  wave  boundary  conditions 
we  obtain  the  eigenfunction  and  wave  eigenfrequency. 
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(ip)  = exp  J-  a(p  -<P0)2/2  I Wn\a^^{\l>  -ipQ 


^ • /)  "a  kc  i 

ci  /d2k  2 s 


o = i -f'  <B^> 
s 


u j q 


1+d 


k%2  ,^'n  ai 


1 ' ^ J ■ <B  > 

Lu  (jJ  N I ' 

e ci 


(2n  + l) 


(15) 


where  are  Hermite  polynomials.  For  d<<l  we  may  approximate  the  shear 

induced  damping  rate  for  the  least  damped  mode  (n=0)  as 

1/2 

Ysh  = - k2cs2  ^ci  1 <B±>  M ' (^)/q(»P)  (16) 

which  will  be  evaluated  in  the  next  section  for  a specific  equilibrium 
configuration. 

We  now  discuss  the  conditions  under  which  the  diagonal  approximation 
to  (10)  is  valid. 

First  of  all,  if  the  plasma  has  a circular  cross  section,  a,n 
and  are  all  independent  of  x , so  that  the  matrix  is  diagonal  because 
of  the  orthogonal i ty  relation.  For  noncircular  cross  section,  the  matrix 
elements  clearly  approach  zero  as  m-p  approaches  infinity.  I he  reason 
is  that  the  matrix  element  is  the  integral  of  a smoothely  varying  function 
(the  poloidal  field)  times  a very  rapidly  oscillating  function. 

In  order  to  truncate  the  matrix  at  a one  by  one,  we  must  invoke  ion 
Landau  damping,  which  is  not  explicitly  included  in  Eq.  (8).  As  the  wave 
propagates  away  from  the  mode  rational  surface,  the  value  of  k ..  increases. 
When  the  distance  from  the  rational  surface  is  large  enough  that  ^/k„v. 


approaches  about  2 or  3,  the  ion  Landau  damping  gets  very  strong.  The 
energy  propagating  away  from  the  mode  rational  surface  is  absorbed  in 
this  region  and  the  eigenfunction  decays  sharply  in  space.  If  the  spatial 
point  at  which  w/k„v.  ~ 3 is  less  than  halfway  between  the  neighboring 
rational  surfaces,  solutions  to  the  diagonal  matrix  (including  ion 
Landau  damping)  are  spatially  isolated  from  one  another.  Thus  the 
eigenfunctions  of  the  diagonalized  matrix  are  a good  approximation  to 
the  solution  of  Eq.  (10). 

The  distance  between  rational  surfaces,  dip,  is 

6^  = (ndq/difO"1  . 


The  distance  at  which  ^/k,, ' 3,&ip,  is  given  by 


A'j;  = 


fLqR__. 


3nv.jdq/d^  . 

where  we  have  used  k=n/R.  Thus  the  condition  that  the  matrix  can  be 
diagonal  reduces  to  <S 2 A . or 


3v  . 

i 

O 

The  work  of  Taylor  applies  to  the  limit  opposite  to  (17) 


(17) 


IV.  EVALUATION  OF  SHEAR  STABILIZATION  FOR  A SPECIFIC  EQUILIBRIUM 

We  consider  an  equilibrium  which  we  partially  specify  by  choosing 
the  axial  current  density  as  a function  of  ip , 


Jz(*>  = Jo  for  °-  ^0  ’ 

(18a) 

Jz(<^)  = 0 for  Ip>i l>0  . 

(18b) 

Thus  Jz  is  a constant  inside  the  flux  surface  ip  = and  is  zero 
outside  this  flux  surface.  Note  that  the  plasma  will  be  assumed  to 
extend  into  the  current  free  region,  > ip  . We  take  the  surface 

ip  = ipo  to  be  an  ellipse  of  ellipticity  k so  that  the  following 
equilibrium  solution  appies  i n ip  < ipQ  (cf.  Ref.  12), 

, kB  2 2 2 2,2  no,\ 

& = 2qR  p » P = x + y /K  » (19a) 

1/2  -1 

= | Vip  | = kB  (x^+y^K  (qR)  (19b) 

This  equilibrium  has  all  magnetic  surfaces  elliptical  with  the  same 
ellipticity.  Also, since  Jz  is  constant  in  ij<  < ti^q  (i|;)  is  also  constant 
in  ^ < tpQ , and  there  is  no  shear.  Outside  the  t[i  = ^ surface  cj ( ip ) is  not 
constant, and  there  is  shear.  We  evaluate  the  shear  induced  damping  in 
ij;  > for  a point  close  to  the  i]i  = tj^  surface.  The  geometry  dependent 
factor  in  (16)  is 


The  evaluation  is  performed  in  Appendix  B for  the  particular  equilibrium 
under  consideration  and  we  obtain  from  (16) 


kc 


-2  il 


sh 


- -2~—r  f ( k ) » f(<)  = i (k  + 1) 


'to  . P 

ci 


3/2 


(21) 


The  variation  of  Y$h  with  <[as  given  by  the  function  f(<)]is  very  weak, 

as  can  be  seen  by  the  following  numbers:  f(l)  = 1 (circular),  f (1.5)  = 0.988, 

f (2. 5)  = 1.104,  f (3)  = 1.242,  f(3.5)  = 1.392.  Thus  for  k between  1 and  2.5 

the  value  of  y h is  within  10%  of  its  value  for  the  circular  case. 

In  particular,  from  a local  theory  calculation  of  the  dissipative  trapped 

13 

electron  mode  stability  in  vertically  elongated  tokamaks  it  was  found 
that  the  growth  rate  was  reduced  by  elongation.  Since  the  shear  damping 
appears  to  be  approximately  independent  of  elongation  (for  this  example), 
this  indicates  the  possibility  of  complete  stabilization  by  vertical 
elongation. 

This  work  was  supported  by  the  United  States  Energy  Research 
and  Development  Administration. 
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APPENDIX  A 


DERIVATION  OF  EQ.  (8) 

Equation  (8)  follows  simply  from  Eqs.  (2)  - (7).  Equations  (2) 
and  (3)  express  ni  in  terms  of  0,  while  Eqs.  (6)  and  (7)  express  vi 
in  terms  of  <T.  Inserting  these  expressions  in  the  ion  continuity 
equation  then  gives  a single  equation  for  4>.  The  only  problem  is  to 
express  gradients,  divergences,  and  unit  vectors  in  flue  co-ordinates. 
The  length  element  is 

ds  = dti,  + dx  + dz.  (Al) 

The  unit  vecotrs  perpendicular  and  parallel  to  the  magnetic  field  are 
then 


ll  = 4 


Bl±z 


(bb  + bb) 


B,i , 
1/2 


i = i 

2r!"  (BB  + Bz2)1/2 


The  components  of  the  gradient  perpendicular  and  parallel  to  the  magnetic 
field  are 


V . = i , B , Trr  + i 9 

~J.  ~{p  i 3iji  ~2 


(i9  • V) 


■ B J + B-L  3Z  ~ °Bz  BJ-  3X 

1 3<l>  (B2  + b2j1/2 


(A2) 
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1 


(A3) 


aB 


V„  - i n 


1 3X 


* Bz  57 


(Bf  + 


Finally,  the  divergence  of  a vector  A is 

A, 


div  A = aB 


2 ( 3 


X 


3t  aB^ 


1 


+ ± ^ + 

3X  B_l  aB2 


JLA 

3z  hz 


( A4 ) 


By  using  Eqs.  (A2)  - (A4)  and  also  assuming  BZ-»B±  (consistent  with 
tokamak  ordering),  and  also  aB^  ^ »kz  (consistent  with  propagation 
nearly  perpendicular  to  B),  Eqs. (2)  through  (7)  can  be  manipulated 
into  Eq.  (8) . 
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APPENDIX  B 


EVALUATION  OF  6$h  FOR  ELLIPTICAL  MODEL 

Since  J =0  in  <p>ip  , VxB  = 0 there  and  consequently  Eq.  (1  ) 

z u * 

implies  that  a is  constant  in  4>  > ip0.  Thus 

d/d^aB^dx]  = -2  f B^3  (dB±/dll(,)dlx 

Also  from  VxB^  = 0 it  follows  that  B^dB^/dl^  = -Rc ~3, 

where  Rc  is  the  radius  of  curvature  of  the  flux  surface  in  the  constant 
z plane.  Thus 


AP  r ° * ri  \ 1/2  fP/R2R  rl  d X H 
a*/'2  (,•/,)-*  h . on 

W?-**/  Jo' 


. dl 
-Bx  dx 


To  evaluate  this  quantity  on  a surface  of  constant  ip  we  note  that  near 
ip  = ipQ  i n »p  > ipQ  the  flux  surface  is  approximately  elliptical  and  B^  is 
approximately  B^  on  the  ip  = ipQ  surface.  To  evaluate  the  integrals 
we  need  B±,  Rc  and  dlx/dx  as  functions  of  x on  the  magentic  surface 
(p  = constant).  From  Eqs.  (19)  we  obtain 


[1  + (k2-1)  (x/p)2] 


1/2 


1/2 


(B2) 


1/2 


_1  /? 

dlx/dx  = [1  + (dy/dx)2]  = [1  + U2-l)  (X/p)2]  [1  - (x/p)2]  . 

o find  R we  express  i . , the  unit  normal  to  the  magentic  surface, 

c ~ v 

Vp  / 1 Vp  I . Rc  is  then  given  by 
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